% function fastPlot(fast)
%   This plots the fast.data in a FAST structure.
% 
% This provides a plot of:
%   All of the binomial trials as well as the current best function fit.
%   The distribution of p(response) at different values along the X-axis
%   (to determine if we were systematically predicting inappropriate
%   values -- a sign that our assumed functional form is wrong)
%   The psychometric function fit.
% 
% This currently plots in linear-linear space, would be nice to add some
% sort of wise decision rule as to when to use semilogx, semilogy, or
% loglog plots... (EV: 1/29/2007)
%
% copyleft Ed Vul & Don MacLeod, 2007
% contact: evul@mit.edu
% version: FAST v2.5

function fastPlot(fast)
% load constants
fastSettings;

if(isempty(fast.data))
    error('fastPlot error: no fast.data present in supplied fast structure');
    plotfast.data = 0;
end

%% plot logarithmically or linearly.
if((min(fast.data(:,1))>0) && ...
    ((log10(max(fast.data(:,1))) - log10(min(fast.data(:,1)))) > ORDERMAGTHRESH)) % logspace
    logx = 1;
    xs = 10.^linspace(log10(min(fast.data(:,1))), log10(max(fast.data(:,1))), 100);
else
    logx = 0;
    xs = linspace(min(fast.data(:,1)), max(fast.data(:,1)), 100);
end

switch(func2str(fast.func.psych));
    case {'psyWeibull', 'psyNormal2', 'psyHalfNormal', 'psyGumbel'} % all functions defined over log 
        logy = 1;
    case {'psyLogistic', 'psyNormal'} % all functions defined over linear space 
        logy = 0;
    otherwise
        warning(sprintf('Psychometric Function %s not recognized by fastPlot, tweak the fastPlot code if you made your own.', fast.params.psych.funcname));
        psyplot = 1;
end

if(logx && logy)
    fplot1 = str2func('loglog');
    fplot2 = str2func('plot'); % semilogx, but no errorbar equiv.
elseif(logx)
    fplot1 = str2func('semilogx');
    fplot2 = str2func('plot');
elseif(logy)
    fplot1 = str2func('semilogy');
    fplot2 = str2func('plot');% semilogx, but no errorbar equiv.
else
    fplot1 = str2func('plot');
    fplot2 = str2func('plot');
end

%% plot fast.data;
figure();
legendkey = {};
subplot(3,1,1);
if(~isempty(fast.data))
    pos = find(fast.data(:,3) == 1);
    if(~isempty(pos))
        fplot1(fast.data(pos,1), fast.data(pos,2),'g.');
        legendkey{end+1} = 'R=1';
        hold on;
    end
    neg = find(fast.data(:,3) == 0);
    if(~isempty(neg))
        fplot1(fast.data(neg,1), fast.data(neg,2),'r.');
        legendkey{end+1} = 'R=0';
        hold on;
    end
    neg = find(fast.data(:,3) == 0.5);
    if(~isempty(neg))
        fplot1(fast.data(neg,1), fast.data(neg,2),'y.');
        legendkey{end+1} = 'R=0.5';
        hold on;
    end
end

colors = {'k -', 'b -', 'c -', 'm -', };
% paramsets = {'margMean'};%{'maxmarginals', 'latticemax', 'initial', 'bestquadfit', 'bestquadfit'};
whichp = fastPsyScale(DEFAULT_PLOT_Ps, fast.params.nchoice, 0);
for i=[1:length(whichp)]
    prediction = squeeze(fastChooseYp(fast, xs, whichp(i)));
    fplot1(xs, prediction, colors{i}, 'LineWidth', 2);
    hold on;
    legendkey{end+1} = sprintf('P = %2.2f', whichp(i));
end
legend(legendkey, 'Location', 'EastOutside');
% ps = '(';
% for i=[1:fast.params.n-1]
%     ps = sprintf('%s %0.5g,', ps, fast.params.est.(paramsets{1}){i});
% end
% ps = sprintf('%s)', ps);
axis tight;
title(sprintf('fast.data and prediction.'));
xlabel('Function X');
ylabel('Function Y');

%%
subplot(3,1,2);

allx = sort(unique(fast.data(:,1)));
for i=[1:length(allx)]
    idx = find(fast.data(:,1) == allx(i));
    px(i) = mean(fast.data(idx,3));
end
plot(allx, px, 'b* ');
hold on;
if(length(allx) >= length(fast.data(:,1))./5) % lump fast.data;
    maxx = max(fast.data(:,1));
    minx = min(fast.data(:,1));
    grx = maxx-minx;
    tenth = grx./10;

    if(tenth > 0) % if we do not have just one x value
        for i=[1:10]
            gd = (((fast.data(:,1)<(minx+tenth.*i+eps)) + (fast.data(:,1)>=(minx+tenth.*(i-1)))));
            idx = find(gd==2);
            x(i) = mean([minx+tenth.*i+eps, (minx+tenth.*(i-1))]);
            stder(i) = 1./sqrt(length(idx));
            mp(i) = mean(fast.data(idx,3));
        end
    else
        x = minx;
        stder = 1./sqrt(length(fast.data(:,3)));
        mp = mean(fast.data(:,3));
    end
    errorbar(x, mp, stder, 'g*');
end
title('P(response) at Each X');
hold on;
plot([min(allx) max(allx)], [.5 .5], 'k -');
ylim([0 1]);
xlim([min(allx)-eps max(allx)+eps]);
xlabel('Function X');
ylabel('P(response)');

%%
xs = [min(fast.data(:,1)):max(fast.data(:,1))];
yscrit = squeeze(fastCalcYs(fast, fast.data(:,1), -1));

if(logy)
    funcvalues = (fast.data(:,2)./yscrit(:));
    org = logspace(log10(min(funcvalues)), log10(max(funcvalues)), 100);
    rg = org .* yscrit(1);
else
    funcvalues = (fast.data(:,2)-yscrit(:));
    org = linspace(min(funcvalues), max(funcvalues), 100);
    rg = org + yscrit(1);
end
            
pp = fast.func.psych(fast.params.nchoice, ...
                     fast.params.est.('marg').mean{fast.params.n}, ...
                     yscrit(1), rg);

subplot(3,1,3);
fplot2(org, pp, 'b -');
hold on;
fulrg = diff([min(funcvalues), max(funcvalues)]);
fplot2(funcvalues+(rand(length(funcvalues),1)-.5).*.02.*fulrg, min(max(fast.data(:,3)+(rand(length(funcvalues),1)-.5).*.02, 0), 1),'kx');

tenths = fulrg./10;
if(tenths > 0)
    for i=[1:10]
     gd = (((funcvalues<(min(funcvalues)+tenths.*i+eps)) + (funcvalues>=(min(funcvalues)+tenths.*(i-1)))));
     idx = find(gd==2);
     x(i) = mean([min(funcvalues)+tenths.*i+eps, (min(funcvalues)+tenths.*(i-1))]);
     stder(i) = 1./sqrt(length(idx));
     mp(i) = mean(fast.data(idx,3));
    end
else
    x = min(funcvalues);
    stder = 1./sqrt(length(fast.data(:,3)));
    mp = mean(fast.data(:,3));
end
errorbar(x, mp, stder, 'g*');
ylim([0 1]);
xlim([min(funcvalues).*.99 max(funcvalues).*1.01]);
legend({sprintf('%s', func2str(fast.func.psych)), 'Data', 'Lumped data'}, 'Location', 'EastOutside');
title('Psychometric fit');
if(logy)
    xlabel('Relative Y (Y/ycrit)');
else
    xlabel('Relative Y (Y-ycrit)');
end
ylabel('P(response)');
